COMPUTATION OF A COMBINED SPHERICAL-ELASTIC AND 
VISCOUS-HALF-SPACE EARTH MODEL FOR ICE SHEET SIMULATION 

ED BUELERS CRAIG S. LINGLE^, AND JED A. KALLEN-BROWN^ 

Abstract. This report starts by describing the continuum model used by Lingle & Clark 
(1985) to approximate the deformation of the earth under changing ice sheet and ocean 
loads. That source considers a single ice stream, but we apply their underlying model to 
continent-scale ice sheet simulation. Their model combines Farrell's (1972) elastic spherical 
earth with a viscous half-space overlain by an elastic plate lithosphere. The latter half-space 
model is derivable from calculations by Cathles (1975). For the elastic spherical earth we 
use Farrell's tabulated Green's function, as do Lingle & Clark. For the half-space model, 
however, we propose and implement a significantly faster numerical strategy, a spectral 
collocation method (Trefethen 2000) based directly on the Fast Fourier Transform. To 
verify this method we compare to an integral formula for a disc load. To compare earth 
models we build an accumulation history from a growing similarity solution from (Bueler, 
et al. 2005) and and simulate the coupled (ice flow)-(earth deformation) system. In the 
case of simple isostasy the exact solution to this system is known. We demonstrate that the 
magnitudes of numerical errors made in approximating the ice-earth system are significantly 
smaller than pairwise differences between several earth models, namely, simple isostasy, the 
current standard model used in ice sheet simulation (Greve 2001, Hagdorn 2003, Zweck 
& Huybrechts 2005), and the Lingle & Clark model. Therefore further efforts to validate 
different earth models used in ice sheet simulations are, not surprisingly, worthwhile. 



1. Two LINEAR EARTH MODELS AND THEIR GrEEN'S FUNCTIONS 

Lingle & Clark (1985) use as their fundamental tools the Green's functions of two different 
linear earth models. The Green's functions for these models are convolved with the load to 
compute (vertical) displacements of the earth's surface. One finds an elastic displacement 
and a viscous displacement given a current load and a load history, respectively, as 
we will explain. The total displacement is then the sum u = + vX at any time. That is, 
the two linear models are superposed. 

The partial differential equations (PDEs) behind these Green's functions are linear. In 
this report we state these PDEs, which is, in the case of the second model, a nontrivial 
accomplishment (see section 13)) . We then approximately solve these PDEs in a demonstra- 
bly efficient manner. First, however, we describe the two models and their sources in the 
literature. 



Draft February 2, 2008. 

^Dept. of Mathematics and Statistics, Univ. of Alaska, Fairbanks AK 99775-6660. Email ffelb@uaf.edu. 
■^Geophysical Institute, Univ. of Alaska, Fairbanks AK 99775-6660. 



2 



BUELER AND OTHERS 



An elastic, self-gravitating spherical earth. The main equations of this model are 
labeled (27) in (Farreh 1972): 

V • r — V {pgs ■ e^) — pVcf) + gV ■ (ps) + ps = 0, 
V2</) = -47rGV-(ps) 

Here p = p{r) is the density of the earth. (We restrict to only radial dependence because 
the densities used by (Farrell 1972) are for "stratified" earths.) Also, g is the acceleration 
of gravity, r is the full stress tensor, and s is the displacement vector field (i.e. the strain 
field) which we seek. The gravitational potential is described below. Note that s = u is 
the velocity field and that s is just the acceleration. 

The first equation comes from formally linearizing the equation of conservation of momen- 
tum. The second equation is for the gravitational potential. The field cj) is the additional 
gravitational potential "on top of" that caused by the undeformed earth and also additional 
to the potential of any masses outside the earth (including the load). In Farrell's words 
"... is the perturbation in the ambient graviational potential 0i plus the potential of any 
externally applied graviational force field (/)2." 

Actually, equation (27) in (Farrell 1972) has "— w^ps" where we have "/os" because Farrell 
only states the Fourier-transformed-in-time equations. We will only be interested in the 
s = 0, equivalently w = case, however, because we are interested in phenomena on the 
scale of years or centuries, unlike Farrell who was interested in tides. 

Consider point forces or disk loads at the surface. For such loads it is natural to use 
spherical coordinates (r, 0, ip) where the z-axis is along the line between the center of the 
earth and the center of the load. Here 6 is the angle between the position vector and the 
z-axis (i.e. the colatitude) and (p is the longitude. The load is at the "north pole." 

For the resulting stratified problem we seek the components of s which do not vanish, 
namely Sr and in s = Sr(r, 0)er + S6)(r, 6)eg. By symmetry the "toroidal" component of the 
strain, s^p, vanishes everywhere. Also we seek the potential (f) = (p{r, 9). The functions Sr,sg,(j) 
are expanded in spherical harmonics with radially-dependent coefficients; see equation (28) 
in (Farrell 1972). The radially-dependent coefficients, for each degree in the expansion, solve 
a system of ODEs in the radial coordinate r. Using a radially-dependent density for the 
earth these can be solved numerically by standard ODE means. Note that (Farrell 1972) and 
(Lingle & Clark 1985) choose the "Gutenberg-Bullen A" model. For the Green's function 
corresponding to a point load, Farrell has done this using a Runge-Kutta method, and we 
accept and use the tabulated result. 

Following (Lingle & Clark 1985) we are only interested in the vertical displacement of the 
surface of the earth, and therefore the vertical displacement u{6) corresponding to a point 
load is the Green's function we seek. In terms of the spherical harmonics expansion the 
relevant equation is equation (37) in (Farrell 1972). Farrell computes this Green's function 
and reports its values at particular distances in his table A3. Table 1 in (Lingle &: Clark 1985) 
also reports this data. One must be clear on normalization so a plot is in order here. Let 
G^{r) be the vertical displacement caused by a 1 kg mass at the north pole and evaluated at 
a distance r along the surface of the earth. (The coordinate r here has a different meaning 
from the spherical coordinate of the same name. Specifically, the new variable r = a9 if a 
is the radius of the earth and 9 is the spherical coordinate, the radian colatitude.) Figure ^ 
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shows G^{r). There is a 1/r singularity to this elastic Green's function, in contrast to the 
Green's function for the flat, viscous model which follows. 




r (km) r (km) 

Figure 1. Two views of the vertical surface displacement Green's function 
G^{r) for the elastic spherical self-gravitating earth model (Farrell 1972). 
Here r is the distance along the surface of the earth from the point of applica- 
tion of the load. Left: the smooth normalized form rG^{r). Right: the same 
data without normalization, suggesting the actual 1/r singularity. Note log 
scale on the horizontal axes. 

This elastic Green's function is used as described in equation (20) in (Lingle & Clark 1985) 
and as follows. Suppose we seek the vertical displacement = u^{x, y, t), caused by elastic 
deformation of the spherical earth. Suppose the load at time t is given by the function 
^{x,y,t), with units of mass per unit area. Then 

u^ix, y,t) = J J G^i\r - r'|)M/(x', y' , t) dx'dy', (1) 
R 

where we define |r — r'p = (x — x')^ + (y — 2/')^> course, and where R denotes a map-plane 
region containing the load. Clearly the displacement depends on time only through the 
changing load; elastic changes are instantaneous. 

In using we necessarily project the earth's geoid into a fixed plane. This projection 
means our results are limited to an appropriately small region of the earth's surface. We do 
integral numerically as explained in section [21 

A viscous, fiat earth overlain by an elastic plate. Next we describe the time-dependent 
Green's function for a model which comes from (Cathles 1975). The PDE actually solved by 
this Green's function is given in section |21 

Cathles' sub-subsection III. A. 2. e, pp. 50-55, describes a viscous half-space asthenosphere 
overlain by an elastic plate lithosphere. An important point about this model, which partly 
explains the superposition "u^-l-u^" used by (Lingle & Clark 1985), is that the elastic plate 
lithosphere used here deflects but does not compress in the vertical. Therefore all vertical 
motion in this model is really asthenosphere motion, though the elastic plate spreads the 
influence of any load. The just-described spherical elastic earth exhibits elastic compression, 
however. 
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Unfortunately, all we are given in (Cathles 1975) are the Hankel transforms of the actual 
equations. Roughly speaking, from (Cathles 1975) we use equation (III-35) along with the 
definitions of "a" and "D" contained in the footnote on page 52 (Lingle 2005, personal 
communication). Precisely speaking, however, our source is (Lingle & Clark 1985), from 
which we use equations (4), (7), and (8)~(14). We will also use the particular choices of layer 
thickness, viscosity, and flexural rigidity for the "two-layer" model from that source. We now 
repeat some equations from (Lingle & Clark 1985) as needed for clarity. 

Let M^(r, t) be the vertical displacement of the surface supposing a point load at the origin 
r = applied at time t = and held. Consider the Hankel transform of this function 

-.V/ 



u {K,,t) = / u {r,t)Jo{Kr)rdr, 
Jo 

where Jq is the Bessel function of zero order (see Appendix^. The Hankel transform is 
self-inverse, so can be recovered from vX by the same integral. 

The half-space model hypothesizes (Lingle &i Clark 1985) that vX solves the equation 

dt 2r]K ^ ~ 2riK ^ ' 

where 

/ , Dk'^ , ^ ET^ 

ain) = 1 H and D = — kt. 

^ ' Prg 12(1 - ly^) 

We denote by azz the Hankel transform of the normal stress from a point load applied at the 

origin 

azz{K,t) = -^H{t), (3) 

in units of Nm~^, corresponding to a point mass of 1 kg. Here H{t) is the Heaviside function 
{H{t) = 1 for t > and H{t) = otherwise). Note ^zzi^^^t) is the Hankel transform of 
(^zzij^t) = —g6o{r)H(t) where So{r) is the Dirac delta function at the origin which acts on 
functions on the plane. ^ The initial condition to @ is the condition of zero displacement 

n^(K,0) = 0. 

That is, vy{K,t) is the Hankel transform of the Heaviside Green's function of an as yet 
unstated PDE of which Q is the Hankel-transformed version; see section 3 for clarification 
of this description. 

Poisson's ratio and Young's modulus for the elastic plate lithosphere are assumed to be 
u = 0.5 and E = 6.6 x 10^'^ M/m^, respectively. The lithosphere thickness T is assumed to 
be 88 km. The resulting flexural rigidity is D = 5.0 x 10^^ N m. The density and viscosity of 
the fluid in the underlying half-space are assumed to be pr = 3300 kg/m^ and t] = 10^^ Pas, 
respectively. 

Equation is an uncoupled set of linear flrst order ODEs in time. That is, the spatial 
Hankel transform has done its job and turned a PDE into a solvable system. Let'^ = 
Prga(K) = pr-g + Dk^. The solution of © and © is 

2ttI3{k) 

^See Appendix^ and especially equation 14211 for the defining property of So{ 
1985) 



■^Use of f3 instead of a represents an admittedly minor simplification of the notation in (Lingle & Clark 
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for t > and u {K,t) = for t < 0. Because of the self-inverse property of the Hankel 
transform, we have the following integral formula for the Green's function: 

G^{r,t) = u^{r,t) = —^j^ /3(K)-i(^l-exp[-^(K)V(2r?K)]) Jo(rK)KdK, (5) 

for t > and G^{r,t) = for t < 0. This formula is equation (14) in (Lingle & Clark 1985). 
Note has units m kg ^; see formula © below. 

As far as we know the integral © must be computed numerically. Furthermore there 
seems to be no one-dimensional procedure analogous to the Fast Fourier Transform (FFT) 
(Bracewell 1978) to do the job quickly. Our strategy for the similar disc load integral (Ap- 
pendix^ is to break up the oscillatory integral into more than 100 subintervals and call 
an adaptive quadrature routine for each subinterval. This strategy is essentially the same as 
that described on page 1104 of (Lingle & Clark 1985) for Q. 

Unlike the elastic case, the Green's function is time-dependent. A graph of G^ for 
several t values is shown in figure |2l The viscous behavior is clear, as is the role of the 
elastic plate lithosphere in removing any singularity at r = 0. Note that the peripheral bulge 
develops only at large times; compare the "standard" model in section [l] 




r (km) 



Figure 2. Green's function G^{r,t) for the viscous flat earth model. The 
curves are from times given in table 1 in (Lingle &: Clark 1985); compare 
figure figure 5(a) there. The top curve is at t = 20 years and the bottom at 
10^ years. Not log scale on the horizontal axis. 

A method for using G^ to compute the response to arbitrary load is described in equations 
(19) and (22) in (Lingle & Clark 1985), and as follows. Suppose there is a load function 
^(x, y, t), with units of mass per unit area, on some region R of our map-plane. The time rate 
of change of this load function, or, equivalently, the incremental changes in this load function, 
are what we sum using G^ to find the general (viscous) displacement vX = vY {x, y, t) of the 
surface. In fact, let 

^ 

Hx,y,t) = —{x,y,t). 
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Then Lingle & Clark assert (with a heuristic motivation; pp. 1105-1106) that: 

u^{x,y,t)= j ff G^{\r-r'\,t-t')X{x',y',t')dx'dy'dt'. (6) 

J —CO J J 

R 

Note that the right side of © has units of length because has units mkg~^, as noted, 
and A has units kgm^^ s^^. 

Lingle & Clark (1985) describe the function A discretely, thereby incorporating an approx- 
imation of the rate of change of the load. Some such approximation is essential in practice, 
of course. They suppose a fixed sequence of past times {ti} and define 

rU 

A(x, y, ti) = / A(x, y, t) dt = ^'(x, y, U) - ^'(x, y, ti_i), 

Jti-l 

so 

u^{x,y,t)^Y. jj G''{\r-v'\,t-ti)A{x',y',ti)dx'dy'- (7) 

* R 

compare equations (16) and (17) in (Lingle & Clark 1985). 



2. The load response matrix method 

Suppose now that the map-plane region R is divided into a grid of rectangular elements of 
area AxAy. Concretely, suppose R = [—Lx,Lx] x [—Ly,Ly] is a rectangular region, suppose 
Nx,Ny are positive integers, and let Ax = {2Lx)/Nx, Ay = {2Ly)/Ny. Center each element 
at {xj,yk) = {—Lx + (j — l/2)Ax, —Ly + {k — 1/2) Ay). There are M = N^Ny elements, each 
denoted by a pair (j. A;) for 1 < j < A''^., 1 < < A^j^. 

We define the elastic load response matrix (LRM) {^{jk){mn)}j j^i^ = '^^■■■,^x, k,n = 
1, . . . , Ny, as the vertical displacement of element (j, k) caused by a unit change in ice thick- 
ness within element {m,n); compare (Lingle & Clark 1985) page 1106. This displacement is 
assumed constant within element {j,k). (As it stands, this description applies only to the 
elastic spherical model in the previous section. Slightly different conventions apply to the 
LRM for the viscous earth model; see below.) 

We compute {^(jk){mn)} by integrating over element {m,n): 

Jyn-Ay/2 Jxm-Ax/2 \ ^ J 

where pi = 910 kg m~^ is the density of ice. Note that the load per unit area (i.e. ^ in the 
last section) for a unit (i.e. one meter) change in ice thickness on a rectangle of area AxAy 
is 

mass Pi(\ ■ Ax ■ Ay) 
area AxAy 

That is, equation Q is a special case of equation 

The matrix {^{jk){mn)} is M x M if we linearly order all elements; alternatively T could 
be regarded as a "4-tensor" with indices j, k,m,n. There are symmetries in this object, and 
we may exploit them to compute roughly {2Nj.){2Ny) = 4M integrals © rather than doing 
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an integral for each of the entries of T^j!. In fact, if we change variables in Q by 
X = Xm- y = Vn- C then we get 

l-Ay/2 /.Ax/2 

rofc)(„„) =p^ / ( - m)^x - + ((A; - n)Ay - C)^ dC- 

J -Ay/2 J -Ax/2 ^ ' 

Let 

/.Ay/2 /.Ax/2 ^ 

/^(p, g) = / / (V(pAx-e)2 + (gAy-C)2) di dC (9) 

J -Ay/2 J -Ax/2 ^ ^ 

for -7V^ + 1 < p < iVx - 1, -A^y + 1 < g < iVy - 1. Then 

r(jfc)(mn) = -m,k- n). 

We need only compute the {2Nx — l){2Ny — 1) entries of I^. Integral @ can be done 
numerically, getting values for the integrand G^(---) by interpolation between the values 
computed by Farrell. 

Now let be the average value of the ice thickness H{x,y,t) over element {m,n). 

Following Lingle & Clark we call {^^(mn) j^""^^^ the load vector. (Technically it is a thickness 
vector, actually; if the load is actually liquid water one computes the equivalent thickness 
to give values Hf^^^y) Integral (pQ), which gives the elastic displacement from the load, is 
approximated by 

A^^ Ny 

U^{Xj,yk,t) « ^ ^'['{jk)imn)H(mn) = J2 ^ Pi - m, k - Tl) Hf^rnn)- (10) 

m=l n=l m=l rt=l 

Note that as a straightforward matrix- vector product, (|1()|) requires M = N^Ny scalar mul- 
tiplications to compute the elastic displacement in the {j, k) element, and thus multipli- 
cations are required to update at each timestep. 

A comparable LRM approximation applies to the viscous model. Let {At,}^^ be a de- 
creasing sequence of Nt positive times. These values are time intervals before the present 
time; table 3 in (Lingle & Clark 1985) gives a list of 26 values. Let be the displace- 

ment in element (j, k) caused by a change in ice thickness of 1 meter in element (m, n) at 
time Ati before the current time. From 

l-yn+Ay/2 /.x™+Ax/2 / , \ 

^\mmn)=Pi / / G^ J(x,-x)2 + (yfc-y)2,At, dxdy. (11) 

Jy-a-Ay/2 Jxm-Ax/2 V ^ / 

Compare to equation Q for the elastic case. Let A^^^-j be the average value of the change 
in ice thickness H{x., y,t — Atj_|_i) — H{x, y,t — Ati) over element (m, n). Then 

Nt Ny 

U^{Xj,yk,t) -J^^Yl %fc){mn)^(mn)- (12) 
1=1 m=l n=l 

As it stands, matrix-vector product requires NfM = N^N^Ny scalar multiplications to 
compute the {j,k) element, thus NtM"^ multiplications to update vY . 

As noted at the beginning, we superpose the results from these elastic and viscous LRM 
approximations : 

m,n i m,,n 

Compare equation (25) in (Lingle & Clark 1985). 
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We are concerned with computability in reasonable time. By using to update every 
element at a given simulation timestep requires 

{l + Nt)-M-M = 0{NtM^) = 0{NtNlNl) (14) 

scalar multiplications. Lingle & Clark (1985) used M = 38, as this was the number of 50 km 
long elements into which the flowline — a single ice stream and included ocean — was divided. 
For Antarctica simulations with Ax on the order of 50 km or so, the minimum reasonable 
number of elements is = Ny at least 80 so M > 6400. Note that Nt is roughly independent 
of the nature of the problem, as long as it involves large amounts of polar ice. Thus by (|14j) 
an Antarctica problem is roughly a factor of 

(64002)/(382) 3 X 10^ 

times more expensive in computation than the problem addressed in (Lingle & Clark 1985) if 
one directly implements the load response matrix method using matrix multiplication. Such 
poor scaling of this numerical method clearly represents a danger when using it in an ice 
sheet simulation. 

A major speedup is possible if one uses the convolution sum form of the multiplication, 
as in equation (|1()|) when is used. (An obvious corresponding construction of a function 
I^{p,q,t) is needed to make the viscous LRM computations above into convolutions. As it 
turns out, we will not need that construction.) Convolutions sums can be quickly computed 
by the Fast Fourier Transform (FFT). Using standard estimates on the time for the FFT 
(Briggs & Henson 1995), we can reduce the time to 

(1 + Nt)-M log M = 0{NtN,Ny log iV, log Ny). (15) 

Compare this with H14|) above. One must still, however, precompute the LRMs, which turns 
out to be more expensive than solving the whole problem (over quite long time scales) if one 
uses the method of the next section. In addition, the method of the next section eliminates 
a factor of Nt work in a time-dependent simulation. The method of the next section also 
significantly reduces memory usage. 

3. The straight-from-the-PDE method 

Derivation. We now reverse engineer some of the Green's function and Hankel transform 
"thinking" in the previous sections. We recover the PDF underlying the half-space viscous 
model (HJ. Actually, the resulting equation is not, technically, a partial differential equation. 
It is a linear pseudo-differential equation easily understood through the Fourier transform. 
We have already computed solutions of this PDF by the Hankel transform, by the indirect 
method of Green's functions. In any case, analyzing the new PDF will lead to a much more 
efficient method for computing deformation in the half-space model. We must still use, for 
now, the Green's function and LRM for the elastic response computed from the Farrell (1972) 
spherical earth model; we will implement the convolution sum ()10p by the FFT. 

Returning to equation (jj}, we apply the inverse Hankel transform. In fact Q is equivalent 

to 

d 

— {2r]Ku) + prgu + DK^u = azz-, (16) 

denoting u = vY for the rest of this section, and with the Hankel transform u = vX . As 
shown in Appendix^ the multiplication by and which appear in equation ()16() can be 
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regarded as the action of operators which are powers of the Laplacian operator. In particular, 

A - -V2 - - (— —) 
acts on the Hankel transform of a radial function f = f{r) by multiplication by k^: 

Af{K) = K^f{K). 

In fact A is the positive Laplacian as we see it is equivalent to multiplication by a nonnegative 
factor. The inverse Hankel transform of (|16j) is 

— (27jA^/^uj+Prgu + DA\ = a,, (17) 

for u{r,t). Equation H17() is the "underlying PDE" for equation (j2J. The symbol A^ stands 
for the standard biharmonic fourth-order differential operator 

^ f — fxxxx ~l~ '^fxxyy ~l~ fyyyy 

(Sneddon 1951, section 20). The operator A^/^ is not a differential operator but is definable 
via the Fourier transform in general; see Appendix \K\ One can also write equation ()17() as 

d 

— {2ri I V| u) + prgu + DV^u = a^^ 
ot 

if the meaning |V| = V— is understood. 

To confirm the equivalence of © and H17|) the reader may verify that the Green's function 
G^{r,t) defined by © satisfies 

^ (27?A1/2G^) +^,5^^ +I)A2g^ = -g6o{r)H{t). 

From now on we remove the assumption of radial load, and suppose equation (|17() applies 
for any load o"^^(x, y, t). The solution n(x, y, t) is a function of three variables; it is no longer 
radial. 

Note that the equilibrium of (|17() is a standard rigid plate equation with a bouyant restor- 
ing force: 

DV^u = 0-22 - PrQU. (18) 

For example, this is equation (8.7.3) in (van der Veen 1999). 

The interesting part of (|17() is the time-derivative term. This term accounts for viscous 
flow within the mantle. It is not completely clear to the author at the present time why the 
particular power A^/^ appears or why it represents the correct diffusive behavior. The units 
in equation (|17() are consistent only with the 1/2 power of the Laplacian, however. 

The Fourier transform of 1)171) is worth noting, because, as in the special case of a radial 
load wherein we may use the Hankel transform, we can write the solution as an integral. 
Namely, if u{^, t) is the two- variable, spatial Fourier transform of u{x, y, t), 

C,t) = - / n{x, y, t)e-*(-'«+^^) dx dy, 

then (|T7)) is equivalent to 

I (2r? {e + ^) + Pr9n + D{e + e)'u = a,,. (19) 
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As with equation ((21), this is a decoupled system of first order ODEs in time. Using initial 
condition u{x,y,to) the solution is 

c,t)= I 2r?(g2+c2)i/2 ^) 



to 

+ exp 



where /3(^,C) = Prd + + C^)^- To find u itself one needs to do the inverse Fourier 
transform, and this could, potentially, be done by the FFT. In fact we will compute more 
directly with H17|) . and we will avoid the integral over time in ()2U() . 

Appendix [0 illustrates the relationship of formula (|2nj) to the Hankel transform formulas 
in section n 

Implementation. We now treat PDF (|17|) numerically by discretizing using a finite differ- 
ence method in time and then computing the action of A^/^ and using the FFT. Our 
method produces a "whole new ball game" numerically relative to integral formulations (the 
LRM method). The resulting new method can be called a Fourier spectral collocation method 
(Trefethen 2000). 

We discretize in time by the trapezoid rule — analogous to the Crank-Nicolson method 
for the heat equation (Morton &; Mayers 1994) — and get an unconditionally stable O(At^) 
method for equation H17|) . In particular, let i„ = nAt for n = 0, 1, 2, 3, . . . and let U^''{x, y) 
be our approximation of u{x,y,tn)- Equation (|17j) is approximated by 

2r] Ai/2 ijn+i^ ^ ^(p^^f/n+i + £)a2c/"+1) (21) 

2r, Ai/2 - ^{prgU^ + DA^U^ + At a,,(x, y,t*). 

Here either azzix,y,t*) = azz{x,y,{n + \/2)At) if the load is known at the time t* = 
{n + l/2)At or azz{x,y,t*) = i {<7zz{x,y,tn) + ^^^(x, y, t„+i)) if the load is only known at 
the times t„, both choices preserve O(At^) accuracy and unconditional stability. 

Equation H17() needs boundary conditions and, in fact, we assume u{x, y, t) 0, and 
similarly for a sufficient number of its derivatives, as (x, y) oo. We also assume that the 
support of the continuous function azz{x,y,t) is bounded for each t. That is, we assume 
there is a zero boundary condition at infinity for the rigid plate and that the load is zero at 
sufficient distance from the origin. 

With careful attention to the boundary condition at infinity, our PDE in its time-discretized 
form, namely equation (|21j) . can be well-approximated by its discrete Fourier transform 
(DFT) version.^ A very reasonable way to incorporate the DFT is to assume periodicity 
in the spatial variables.^ For convenience we will also assume a square region. In fact, we 
assume L is the half-length of a computational domain (x, y) £ Q = [—L, L] x [—L, L]. This 
domain may be substantially larger in practice than the desired region of physical interest 
[—Lx, Lx] X [—Ly, Ly]. We will apply periodic boundary conditions at x, y = itL and we want 
L to act like oo when we do this. See figure |31 



"^The discrete Fourier transform is the name of the mathematical operation; the FFT is an algorithm for 
computing the DFT (Briggs & Henson 1995). 

^Other boundary conditions could be applied along the boundary of 57 — e.g. a clamped condition — but 
none of the easily implementable choices are obviously superior. 
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Figure 3. Fast Fourier Transform methods require periodicity. We impose 
periodic boundary conditions significantly far outside the region of modelled 
load changes or significant deformation. 



The PDF problem to which we apply the DFT is, therefore, Fquation (|21|) on the interior 
of O = [—L, L] X [—L, L] with periodic boundary conditions on the boundary of and with 
the initial condition that U^{x,y) is known on Q. 

Note that from (|21|) . if L = oo then the (non-discrete) Fourier transform [/" = J^2U^ 
satisfies iteration 

fjn+i . , .^ ^ [2r?Ac - (At/2) {prg + Dk^)] C) + At a^^, C, t*) 

2rjK + {At/2){prg + DK^) ^ ' 

where = + C^- This is an easily-computed iteration if one can do the (non-discrete) 
Fourier transform T2 exactly. 

To use the DFT we transform the problem to a standard region Cl = [— vr, vr] x [— vr, vr]. Let 
X = -Kx/L, Y = ny/L. Then (|2T|) is equivalent to 

(2r//x + ^(p^5[/"+i + Z)/a2;7"+1) (23) 

= (2w.aV2[/-) _^(p^g[/n + D/A2;7") + At^7,,(X,y,t*) 

on n, where = C/"(X,y), fi = vr/L, and A = - [d'^/dX'^ + d^/dV^). 

Let N be an integer; typically is a power of 2 for efficiency in the FFT. Let h = 271 /N 
and let Xj = —ir + jh, = — vr -|- kh, j,k = 1, . . . ,N. On we use the DFT, as normalized 
by Trefethen (2000), in variables X,Y. If f{X,Y) is some function on 0, with grid values 
fjk = fi^jj^k) then the DFT (forward, inverse) pair is 

N N/2 

fp, = h'J2 e-^^^^^+'^^^V,., fjk = E ^'^'''^^''^'^fpr (24) 

j,k=l ^ ' p,q=~N/2+l 
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Let 

N/2 

fiX,Y) = --^ e^(^^+^^)/p. (25) 

^ ' p,q=-N/2+l 

be the "band-limited trigonometric interpolant" of f{X,Y) (Trefethen 2000); note the rela- 
tion to the inverse DFT (PH) . We see that 

N/2 

^ ' p,q=-N/2+l 

That is, the Laplacian A on corresponds to multiplying the p, q mode by (p'^ + q'^)- Thereby 
A^/^ and A^ are also defined, respectively, by multiplication by (p^ -|- g^)^/^ and (p^ -|- q^)^. 

It now follows that (|23() is very easy to compute if we approximate U^{X, Y) by its band- 
limited interpolant JJ'^{X,Y) and compute the action of the powers of the Laplacian by 
the multipliers above. This describes a Fourier spectral collocation method. That is, one 
time-step in solving PDE (|17() by our method is the sequence 

(i) compute the DFT U^g by FFT from values ^ U''{Xj,Yk) = U''{xj,yk); also 
compute the DFT of the load {azz)pq at t = t* from values crzz{xj,yk,t*), 

(ii) compute 

^ [2r?/x(p^ + g^) - (At/2) (p,g + D^,\p' + q'f)] U^^ + At (a,,)^^ 
""^ 27]fi{p^+q^y/^ + {At/2){prg + D^i^{p^+q'^y) ' ^ ' 

where p, = tt/L and 

(iii) undo the DFT (i.e. do the inverse FFT and make sure the result is real) to get 
Compare equation (|26|) to equation (|22|) which applies for the non-discrete Fourier transform 

Full Matlab implementations of the methods in this report are given in Appendix ^ 
Only a few lines of Matlab are needed to implement the core sequence above, however: 

for n=0:M-l 

[computations using current displacement uun = U'^{xj,yk)] 
uun=uunl ; 

[get E = H{xj,yk,t*)] 
sszz=-rhoi*g*H; 

frhs=right.*fft2(uun) + f f t2(dt*sszz) ; 
uunl=real(ifft2( frhs./left )); 
end 

Here "right" and "left" are pre-computed grid values of the expressions 2rjfik — {At/2)B 
and 27] fik + {At/2)B which appear in (|26|1 . respectively, where k'^ = p^ + q^ and B = 

PrQ + Dfi'^k^. 

We call the iteration (i), (ii), (Hi) the "PDE method" in contrast to the load response 
matrix method ("LRM method") of the previous section. With standard estimates on the 
speed of the FFT when A^^^, Ny are powers of two (Briggs & Henson 1995), the "PDE method" 
requires 

0{N,Ny{log2N,){log2Ny)) (27) 
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scalar operations to update the vertical displacement. This compares directly to equations 
((TH) and (dSI) for the "LRM method." In particular, relative to the 0{NtNxNy{log2 iVa;)(log2 Ny)) 
estimate (jl5|) for the "LRM method" using the FFT for convolution we note a factor of Nt 
less work. Furthermore, the "PDF method" avoids the entire stage of computing the LRM 
integrals, which turn out to be quite expensive, though totally uninteresting, computations. 

Relative to the "LRM method" without the FFT, using representative values of Nt = 25 
and Nx = Ny = 80, and assuming that the constants in the "big O notation" are about the 
same, we get a speed up of about 

(25 • 80^ • 80^)/(802(log2 80)2) « 4 x 10^ 

This ratio is roughly what we observe in practice. For instance, on the same computer 
we compared Fortran 77 codes running the "LRM method" sans FFT using Nt = 101 and 
Nx = Ny = 31 with Matlab codes (Appendix^ implementing the "PDF method" using 
A^^ = Ny = 32. We used At = 500 years and bed deformations were computed for 50k years 
in both cases. The former method took about 10 hours while the later took about 4 seconds 
for an observed speedup of about 9 x 10^. 



4. Discussion 

In the next section we describe the results of computations with the "PDF method." It is 
appropriate, however, to first directly address the apparently new feature in PDF ()17|). 

^ (27] n) + prgu + DA'^u = a,^, 

namely the viscosity expression ''d/dt (2r?Ai/2u)." 

Within the ice sheet modeling community there is a simplified existing standard model for 
aflat "elastic plate (lithosphere) that overlies a viscous asthenosphere" (Greve 2001, Hagdorn 
2003, Zweck &: Huybrechts 2005). Comparison to this model illuminates the significance of 
the viscosity expression in H17|) . The standard model consists of two equations 

Prgw + DA'^w = (28) 

-dt - 7 ' ^^^> 

where u^{x,y,t) is the vertical displacement of the earth's surface, UQ{x,y) is a hypothesized 
unloaded displacement state, and w{x,y,t) is the position of a notional elastic plate which 
is in equilibrium with the current load azzix,y,t). The essential viscous constant for the 
standard model is a characteristic time scale r of relaxation, chosen, for example, as 3000 
years by Greve (2001) and Zweck & Huybrechts (2005). The relaxation time r is indirectly 
related to the asthenosphere viscosity rj; more on this below. 

We can now calculate an illuminating comparison by hand. Suppose that at time t = all 
load is removed but that the vertical displacement is a y-independent sinusoidal mode with 
spatial frequency k: 

u{x, y,t = 0) = u^{x, y,t = 0) = Aq ex.p{ikTTx/ L). (30) 



Here is the initial amplitude and L is a characteristic length scale. We ask: how does 
such a mode decay in the two models? 
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Note that the qth power of the Laplacian act on this mode as fohows: 
A^expiikirx/L) = {k^n'^ / L'^)" exp{ikTTx/ L). 
Thus model (jlZj) has solution u{x,y,t) = A{t) exp{ikTTx/L) where 

2r]{\k\Tr/L)A + prgA + D{k^7r^/L^)A = 0. 
That is, in model H17|) the amplitude of the k mode satisfies 

■ prgL^ + Dk'^Tr'^ . . . ^ . 

^ = - 2,L>m. ^' = pi) 

so a mode with frequency k decays a rate that depends upon k. By contrast, in model (|28|) . 
H29|) with no = the same mode evolves by 

A = --A, A{t = Q) = Ao (32) 
r 

because w = as the load has been removed; equation reduces to du^/dt = —t~^u^. 
Here we see that all modes decay at a rate independent of the frequency. 

But it is clearly the case that a viscous asthenosphere will make elastic plate modes decay 
at different rates depending on the frequency. Indeed, Greve (1997) identifies the failure of 
the standard model H28|) . (|29() to have frequency dependent relaxation times as a deficiency 
of that model relative to full spherical self-gravitating models. Comparing equation (|31|) to 
we are motivated to plot the function 

2riL^\k\-K , , 

which has units of time. Supposing L = 2000 km and that pr,g,D,rj have the values given 
in section^ we plot r(/c) in figure 0] We see that the standard choice r = 3000 a in 
corresponds to frequencies A: ~ 1 and A; ~ 10 in (|31j) . but that no constant relaxation time is 
representative of the actual relaxation spectrum. 

Furthermore we see in formula (|33|) and figure|l]that for small k (e.g. A; < 3) the relaxation 
time T{k) is proportional to k. This behavior is identified by (Klemann & Wolf 1999) — 
see Greve (2001, section 5) — as correct for the most significant mode in a spherical, self- 
gravitating viscoelastic earth model. 

The justification for model (|^ . is its computability, of course. Indeed, the com- 
putation of the elliptic PDE ()28|) is standard in all numerical paradigms (finite difference, 
finite element, spectral). At least on a rectangular spatial grid, however, the time-semi- 
discretization (|^T|) of equation (fTTj) is just as computable as (j2Sl)- In particular, if an ice 
sheet simulation is performed on a rectangular grid using a finite difference or finite ele- 
ment method for the ice dynamics then equation (|17|) can be easily computed by the Fourier 
spectral collocation method of the previous section. 

Now we come to another reason to prefer equation H17() as a model for earth deformation 
in the context of ice sheet modeling. Let us suppose that at the current time the ice thickness 
Hq (and possibly water depth, giving an effective ice thickness) in a region of interest has 
been well-measured. Let us suppose that a reasonably detailed map of current uplift rate 
iiQ = du/dt is also known. This is a realistic supposition given given the state of observational 
geophysics circa 2006 because uplift can be well-constrained by GPS measurements (Larsen 
and others 2005) when bedrock is exposed. Alternatively a spherical viscoelastic earth model 



COMPUTATION OF A VISCOELASTIC DEFORMABLE EARTH MODEL 



15 



11000 




frequency k 

Figure 4. Frequency dependent relaxation time T{k) (solid) for mode k in 
equation (|3Ujl . The value r = 3000 a (dashed) for the standard model is a 
reasonable constant value, but no constant provides a good fit. 

of more-or-less arbitrary sophistication and computational expense might generate a trusted 
current uplift map (Ivins &: James 2005). In either case we can then use ()17|) to determine the 
initial condition for the earth deformation part of an ice sheet simulation without requiring 
further reference to an assumed past load history; compare the integration over load history 
scheme used in (Lingle & Clark 1985). In fact, by p7j) we may solve 

Prguo + DA\o = PigHo - 2riA^/^uo (34) 

for no to get the starting displacement. In other words we ask for the "pre-bent" position of 
the elastic plate in the half space model which accounts for the current uplift rate using the 
current load (i.e. current effective ice thickness). Solving H34() numerically is no harder than, 
in fact it amounts to, one step of the numerical method already described. If the load does 
not change in the simulation, an uninteresting case for ice sheets of course, then the elastic 
plate overlying the viscous half space will start at the current time with the current uplift 
but will then relax to the state satisfying equation (|18|) for equilibrium with the load and 
bouyant force, and there will no more uplift. Note that the presence of bed topography is 
completely irrelevant here because of the linear nature of the model; see comments in (Lingle 
& Clark 1985) to the effect that bed topography represents a irregular "thin veneer of zero 
strength" atop an elastic plate lithosphere of (significant) flexural rigidity D. 

The mechanism described in the previous paragraph is, we believe, a more principled 
replacement for the hypothesized "unloaded surface elevation" UQ{x,y) used in the standard 
model ((^ . Use of that standard model seems to require an assumption of present 

day isostatic equilibrium with the current load (Zweck & Huybrechts 2005) or other artificial 
assumption which is in conflict with observed spatially- varying current uplift. 

An entirely different class of viscoelastic earth models exists in the literature, of course. 
These are the layered, spherical, self-graviating viscoelastic models which typically descend 
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from the work of Cathles (1975) and/or Peltier (1974). The numerical implementation of 
these models typically involves computing a high degree spherical harmonic expansion of the 
strain field for the entire three-dimensional geoid. The traditional difficulty with these models 
is their computational cost (Fastook 1999, Greve 2001). Furthermore there is only modest 
benefit because results from the standard model ((^ . (f^^ above, in particular, are regarded 
by the ice sheet modeling community as reasonably close to those from the spherical models 
(Greve 2001). Equation (|17j) is promising because it is just as computationally inexpensive as 
the standard model but incorporates at least one important feature of the spherical models, 
namely frequency dependent relaxation times. 

Speaking mathematically, an interesting additional possibility exists within the same class 
of computationally inexpensive "PDF methods." Namely, one should be able to modify 
equation (|17j) to take spherical effects into account. In particular, equation (|17|) would be 
computable at essentially the same speed if it were replaced by a non-constant coefficient 
version, for instance 

^ (2r]{x, y)Ai/2^) + prgaix, y)u + D{x, y)A {l3{x, y)Au) = a,,. (35) 

We do not currently know that a set of non-constant coefficients r],D,a,(5 exist which cor- 
rectly account for spherical geometry. It seems, however, that classical continuum mechanics 
and differential geometry must produce such a form because the standard model for a spheri- 
cal self-gravitating earth is a linear model. That is, it has linear response to load. Abstractly, 
this linearity is all that is necessary to make it inevitable^ that a two (spatial) dimension, non- 
constant coefficient linear equation for the vertical displacement of the earth's surface must 
exist for each patch of the earth's surface. It may well involve additional pseudo-differential 
operators not present in the putative form (|35j) . however. 

5. Numerical results 

A "tweak" to the procedure. It turns out that a small error can be avoided if the "PDF 
method" is modified slightly. In fact, in verifying the "PDF method" below, using an exact 
integral formula for a disc load, we observed that there was a uniform a error of several 
meters. This uniform error decayed slowly as the distance at which the periodic boundary 
condition was applied went to infinity. 

The following simple modification eliminates this error. Using the notation of section El 
consider U^{xj,yk) on the grid at timestep n. Let be the average value of C/" along the 
boundary of the computational domain 17 = [—L,L] x [—L,L]; recall 0, is typically larger 
than the physical region [—Lx,Lx] x [—Ly,Ly]. Let u'^^ Ro(^) t>e the vertical displacement 
at distance r from the center of an ice disc load of thickness Hq and radius Rq of an elastic 
plate in equilibrium with the bouyant restoring force. That is, let ^ (r) be the value 
from formula (|44|) in AppendixEl Choose the values Hq, Rq so that the volume ttR^Hq of the 
disc load matches the current (timestep n) load volume, or rather its ice equivalent volume 
if appropriate. 

Our "tweak" replaces the solution U^{xj,yk) at each timestep n with values to which a 
constant shift has been applied: 

[/"'°(x„yfe) = U"ix„yk) - U2 + ^Ho,RoiL). (36) 

^See the discussion of distributions in (Reed & Simon 1980), for instance. 
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That is, we want the "far-field value" produced by the original "PDE method" to be thrown 
out and replaced by the equilibrium plate value with an equivalent disc load. Though the 
volume of the equivalent disc is determined by the current load, one obviously has some free- 
dom in choosing its thickness and radius. We presume that an effort is made to approximate 
the aspect ratio of the actual load, but close matching is not essential. In fact the Green's 
function value would work reasonably well, too. 

Verification. The first concern regarding computations with our earth deformation model 
is verification. In particular, we want to know if numerical results from the "PDE method" , 
with the just-mentioned "tweak," are close to highly-accurate solutions of the continuum 
equation ()17|) . In seeking such solutions we inevitably come to disc loads. Appendix IbI 
addresses this case by the Hankel transform. It yields equation (|43j) . an integral formula for 
the time-dependent radially-symmetric deflection (r, t) resulting from the application (at 
time zero) of a disc load; see figure 1151 The integral must be computed numerically, but 
numerical quadrature is an approximation completely independent of the cartesian-grid- and 
FFT-based "PDE method." We use integral formula (|43() as an "exact" solution, believing 
the accuracy of our numerical integration of H43|) to exceed that of the "PDE method" for 
any achievable grid. 

Our Matlab implementation of the "PDE method" is the function fastearth.m listed 
in AppendixEl Also listed are implementations of the LRM method for the spherical elastic 
earth model (gef orconv.m) and the numerical integration of equation (|43j) (viscdisc .m). 

Let us define a particular numerical experiment. Suppose we use the parameters specified 
in section^ D = 5.0 x 10^^ N m, pr = 3300 kgm~^, and r] = 10^^ Pas. Suppose the disc of ice 
(density pi = 910 kgm^'^) has radius 1000 km and thickness 1000 m. We seek the deflection 
on a square region R, centered on the disc load, of side length 4000 km [Lx = Ly = 2000 
km). Suppose that at time zero the deflection is identically zero, and suppose that the load 
is applied at that time. We calculate the deflection at t = 20k years as computed by the 
"PDE method" and by formula 

There are three numerical parameters of importance for the "PDE method" : 

• N, the number of grid points in each direction, 

• Z, the factor by which the computational domain = [—L, L] x [—L, L] is larger 
than the physical domain [—Lx,Lx] x [—Ly,Ly]; here Lx = Ly, and 

• At, the time step used in approximating the time derivative which occurs in equation 

m- 

Our veriflcation involves showing that as these parameters go to their continuum limits 
(A^ — > oo, Z — > oo. At — > 0) we approach the "exact" solution 

For verification we first fixed At = 100 years and considered the effect of grid refinement 
and of changes to the distance at which the periodic boundary condition was applied. Re- 
garding grid refinement we considered N ranging over powers of two from 2^ = 8 to 2^ = 256. 
Regarding the distance to the periodic boundary condition we imposed the periodic bound- 
ary condition at = L, where L = Z Lx = Z Ly, and we used Z = 1,2,4,8, but we 
quickly discovered that with the above-mentioned "tweak" any value Z > 2 works fine; not 
shown. 

The result for maximum error under grid refinement is shown in figure |SJ This maximum 
error is not the only reasonable measure. As shown in figure for fine grids errors greater 
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than one meter are highly locahzed to certain points just at the edge of the disc load. Note 
realistic ice loads do not have margins as sharp as this disc load. For these reasons among 
others it is reasonable to consider average errors, and we see in figure [3 that the average 
errors are less than 20 cm for = 256, or roughly 0.07% of the compensation depth. 




2 1 I I I I I I 

8 16 32 64 128 256 

N 

Figure 5 . Maximum error made by the "PDF method" relative to the Hankel 
transform integral (|43() . Grid refinement (increasing A^) reduces the max error 
to below 3 m when N = 256. Here Z = 2 and At = 100 a. 
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Figure 6. Spatial distribution of the error when = 256, Z = 2, and 
At = 100 a. Contours of the error |(PDF method) — (equation at 0.5, 

1,2 m. The error is concentrated where the edge of the disc load meets the 
coordinate axes. 
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N 

Figure 7. Same as figure [S] but now average error. Grid refinement (increas- 
ing N) reduces tlie average error to below 20 cm when N = 256. 

Next we compare the effects of spatial grid refinement, increasing A'^, to reduction of time 
stepsize At on the error. See figure |H1 We see that any value of At less than 500 years is 
fine; this is great news for ice sheet simulation. It is not, however, surprising because of the 
relative timescales of ice versus asthenosphere flow. 
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Figure 8. Average error as in figure [7| but with At varying, and for several 
values of N. There is no need for At < 500 a. Spatial grid refinement (in- 
creasing N) is more important to reducing error than is temporal refinement 
(decreasing At). 
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Ice sheet modeling. Earth deformation used in the context of ice sheet modeUng is our 
actual interest. Earth deformation is obviously effected by ice sheet flow — the load moves 
around. Conversely, as the bed deforms the surface slope of the ice sheet changes and this 
effects flow. There is non-trivial coupling. 

The reasonable simplest ice sheet model is the isothermal model with Glen rheology 
(Paterson 1994, Nye 2000). Let h{x, y, t) be the surface elevation of the ice and let H{x, y, t) 
be the ice thickness. The frozen-base isothermal ice sheet equation is the single nonlinear 
diffusive partial differential equation 



where n is the Glen exponent, e.g. n = 3, and F is a constant (typically F = 2{pig)'^ Aq / {n+2) 
where Aq is a softness parameter). If b{x,y,t) is the ice sheet bed elevation — a slight change 
of notation from the rest of the paper — then of course h = b + H . 

As we now show, exact similarity solutions to this equation which incorporate simple 
isostasy (Nye 2000, Bueler et al. 2005; compare Halfar 1983) provide a very nice tool to 
examine coupling to the earth model. They help illuminate the differences among earth 
models. By "simple isostasy" we mean the rule which specifies 



where / is a fixed fraction of the ice thickness (Nye 2000); we will let / = pi/pr = 0.27576 
in our numerical experiments. Since h = b + H , if equation (|38|) applies then /i = (1 — f)H. 

We will compare numerical results for three coupled ice sheet flow/earth deformation 
models: 



Figure IHl shows a result of a coupled simulation for these three models. A detail near the 
margin is shown in figure THii 

In fact, the result shown in figurelHlcame from starting with H = and 6 = at t = and 
using an accumulation history corresponding to the simlarity solution illustrated by figure lTTl 
That is, the accumulation M{x, y, t) comes from equations (9) and (10) in (Bueler et al. 2005), 
using / = 910/3300, X = 5, a = -1, P = 2, Hq = 3600 m, Rq = 750 km, F = 9.0177 x lO'^^ 
m^^ s^^, and with the additional statement M = Bt^^Hx. Note to = 40034 years. In 
addition, at time t = to the accumulation is turned off and so for t > to the exact behavior 
of the solution to Simple is a Halfar-type (Halfar 1983) accumulation-free solution. Thus 
the accumulation history is from a similarity solution to equation (|37|) . incorporating simple 
isostasy, which grows from zero at t = to maximum height at to = 40034 years and spreads 
out from then on, with no loss of volume. 

The importance of such a similarity solution is that it forms an exact continuum solution 
to the Simple model. Therefore we can answer with some precision the question "/low 
do differences resulting from coupling to various earth deformation models compare to the 
numerical errors which occur in ice sheet modeling?" This is an important question. If 
numerical ice sheet errors demonstrably exceed the earth model differences then we should 
be skeptical of any expenditure of effort in the earth modeling direction. Conversely, even 




(37) 



b 



fH 



(38) 



Simple: equations (|37|) and (|5S)) . 

Standard: equations ((HZI), (EHI), (EHl), and b = u"", 

Lingle&Clark: equations (jSZI), (UJ, ^7^, and b = + . 
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Figure 9. Ice sheet on deforming bed, at time 60k years, from three earth 
models Simple, Standard, and Lingle&Clark. View of gridded numerical 
values {Nx = Ny = 192) along the positive x-axis of the grid. 
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Figure 10. See figure IHl detail near the grounded margin. Exact similarity 
solution to the simple isostasy model is added (solid). 

if the differences among earth models are significant, one should report these differences 
relative to the actual magnitude of numerical ice modeling errors. 
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Figure 11. Views of a similarity solution to equation H37() . Left: volume 
over time. Right: margin radius (solid) and dome height (dashed) over time. 
At time to = 40034 a the solution switches from growing (A = 5) to Halfar 
(A = 0), in both cases with simple isostasy (/ = 910/3300). 

Figure ini indeed suggests differences among the coupled models. We ran each model to final 
time t = 60k a. As shown in figure however, all of the models produce the same volume 
at the final time, and indeed at all times; this follows from using the same finite difference 
approach for the ice fiow (as described in (Bueler et al. 2005)) and, of course, the same 
accumulation history M{x, y, t). So the differences can be described by the distributions of ice 
thickness. In figure [T3l we show the maximum and average of the pairwise absolute thickness 
differences |-ffsiMPLE — -^StandardI, etc. (The average differences are over the H > grid 
points under Simple.) We see average thickness differences greater than 10 m between each 
pair. We see that the greatest pairwise difference is between Simple and Lingle&Clark; 
compare figure IHl 

We see a similar picture for bed elevation differences, with Simple versus Standard 
showing somewhat smaller differences, and the comparison Simple versus Lingle&Clark 
again being largest. 

Now, are these differences significant? The answer shown in figure El is yes. With a 
caveat. As noted in (Bueler et al. 2005), ice sheet flow simulations on grids inevitably make 
large thinkness errors near the margin. These errors decay only slowly under grid refinement, 
as can be seen in figure [T5l 

6. Conclusions 

We have seen several modeling and computational issues and numerous equations. So let 
us identify our major point: Equation (fTTjl 



where A is the positive Laplacian A = —d'^/dx^ — d"^ /dy'^, is both 

• a better model for a viscous half space overlain by an elastic plate than the standard 
model (|28|) . (|29|) which is widely used in the ice sheet modeling literature. 
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Figure 12. All models have identical convergence of numerical volume at 
t = 60k a; they share the same accumulation history. 




Figure 13. Maximum and average ice thickness differences in pairwise comparison. 

• and is very computationally tractable on a rectangular grid using a Fourier spectral 
collocation method. 

In brief, one derives equation ()17() by starting with equation (4) in (Lingle & Clark 1985), 
clearing denominators, and then taking the inverse Hankel transform by recognizing powers 
of the positive Laplacian A. That is, equation (fTTj) is equivalent to equation (4) in (Lingle 
& Clark 1985). 
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Figure 14. Maximum and average bed elevation differences in pairwise comparison. 




Figure 15. Ice thickness differences in pairwise comparison as in figure IT^ 
but with numerical errors for the simple isostasy case superimposed. Dif- 
ferences among coupled ice-earth models significantly exceed numerical error 
except for localized numerical errors within a couple of grid points of the 
margin. 
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Of course, equation (|17|) fails to incorporate spherical or self-gravitating effects. Following 
Lingle & Clark (1985) we have, however, chosen to superpose upon the result of ()17() a 
purely-elastic, but spherical and self-gravitating (Farrell 1972), displacement from equation 
This is an admittedly ad hoc way to incorporate spherical and self-gravitating effects 
into an earth deformation model. ^ We note that the effects of equation are "longer 
range" than those of equation (|17j) at large times, so in some sense it is more important to 
incorporate "sphericalness" into the purely elastic part of the earth model. 
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Appendix A. On the Hankel transform and powers of the Laplacian 

A reasonable first reference for the Hankel transform is Chapter 12 of (Bracewell 1978). 
Sneddon (1951) is much more complete, however. Sneddon addresses the application of 
Hankel transforms to elasticity problems, in particular. 

Definition. Suppose /(r) is defined on (0, oo) and suppose that Jg°° \ f{r)\rdr < oo. The 
Hankel transform f{n) , k > 0, is 

l-OO 

f{K) = / f{r)Jo{Kr)rdr. (39) 
Jo 

This is the transform as normalized by (Sneddon 1951), (Cathles 1975), and (Lingle & 
Clark 1985); BraceweU (1978) is shghtly different. 

Here Jq is the zeroth-order Bessel function of first kind 

k=0 ^ ' 



^2k 



an entire function; note | Jo('")| ^ 1 foi" ^ r. J^^z) is the unique solution to the ODE initial 
value problem 

zV(^)+^y'(^) + (^'-0)y(2) =0, 2/(0) = 1, y'(0) = (40) 
and it has integral formula 

1 /■27r 

M^) = :r e-'^'-^'de. (41) 

Jo 

Both (jlOI) and (jH)) will be used below. 

The most natural source of the Hankel transform is as the Fourier transform of a radial 
function on the plane. In particular, suppose /(x, y) is a bounded and integrable function on 
the plane which is actually radial / = f{r). Suppose one computes the two- variable Fourier 
transform J-'2[f] = / by converting the integral to polar coordinates: 



1 



oo ;>oo 



/(?, = ^1 I /(r)e-^(-«+2'^) dxdy= I fir) 



oo 



— oo J — oo 
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1 

2^ 



2lT 



^ — ITK COS 



'dO 



r dr 



fir) 







1 

2^ 



^-ivK cos 9^0 



rdr= / f (r) jQ{nr)r dr 



where = + C^- We have used (|^T|) above. Concretely, we have substituted x = rcosO, 
y = rsinO, ^ = Kcoscj), and ( = Ksin(/). Thus + yC = rKcos{6 — 4>), and for fixed 4> 
the function fi9) = e~*'"'^™^(^~'^) is periodic with period 27r. We conclude that in these 
circumstances / is also radial. 

The Hankel transform (|39|) is evidently linear. The general two- variable Fourier transform, 
which we have normalized to be unitary, has the property 



9ix,y) = g{-x,-y). 
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Therefore the map g g is the identity when restricted to the subspace of radial functions, 
so /(r) = f{r) = f{r) and the Hankel transform is self- inverse. 

For / = f{x,y) sufficiently smooth, define the positive Laplacian of / to be 

We have chosen the sign of the Laplacian to be positive as an operator, which is to say 
that the eigenvalues/spectrum of this operator are nonnegative. In particular, the Fourier 
transform of the Laplacian is multiplication by a nonnegative function 

A7(c,c) = (e' + c')/(c,c). 

For a radial function /(r) we recall that A/(r) = — /" — r~^f'. 
Consider the Hankel transform of the Laplacian. 

Lemma. If f'{r) is bounded as r ^ 0^ and if max{\f{r)\r, \f'{r)\r} — > as r — > +cxd then 
the Hankel transform formula holds: 

Af{K) = K^f{K). 

Proof. Letting prime denote d/dr, two integrations-by-parts do the job: 
( - f"{r) - r'^f'{r))Ja{Kr)rdr 

= f'{r)Jo{Kr)r ^ + / f' (r)— {JoiKr)r) dr - / f'{r)JoiKr) dr 







CO 



oo 



0+ / /'(r) [Jo(Kr)Kr] dr = /(r)Jo(Kr)KrJ^ -y f{r){j'Q{Kr)Kr)' dr 



i - / ^ {z'4(z) + z4{z)) dr^- ^ i-z'Mz)) dr 
Jo ^ Jo ^ 

POO 

= f{r)jQ{Kr)rdr = K^f{K). 

Jo 

We substituted z = nr \n step * to recognize Bessel's equation (jlH]) in step f. □ 

It follows that also for powers of the Laplacian we can give nice Hankel transform formulae. 
In particular, if 

A 2 ^4 d'^ 5^ 

= V = h 2 I 

denotes the biharmonic operator (Sneddon 1951), and if /(r) is radial and has appropriate 
boundedness, then 

Now we define the positive square root A^/^ of the Laplacian A via the Fourier transform:"^ 
Definition. Define (A^/^/)(x,y) by 



'''The Fourier transform gives a spectral resolution of the positive operator A acting on the plane, subject 
to appropriate boundedness at infinity (Reed & Simon 1980). Thus we are using the functional calculus to 
define A^/^ 
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Note that A^/^ acts on scalar functions to produce scalar functions but that it is not a 
true differential operator. In particular, it is neither the gradient nor the divergence. It 
can, however, be easily calculated using the Hankel transform because the following formula 
applies for radial functions / = f{r): 



(Ai/2/)(^) = ^/(^). 

Finally we consider the Hankel transform of the Dirac delta "function." Care must be 
taken because the delta function we need is a function only of the radial coordinate while 
its defining property applies to functions in the plane. In fact, let 6(^xo,yo)i^^y) defined by 
the integral 

/oo fOO 
/ f {x,y)\xo,yo){x,y) dx dy = f{xo,yo) 
-oo J —oo 

for all continuous f{x,y). Let 5o = S{o,o) denote the delta function at the origin. In polar 
coordinates, and for radial functions / = f{r), 5o has the property 

f oo /•2tt 

fir)6oir,0)rdrde = fiO), 



which simplifies to 

/(O) 



f{r)6o{r)rdr = ^-^ (42) 
zvr 



because 5o is independent of 9. Thus 

6o{k) = I 5o{r)Jo{Kr)rdr - '^^^^^ - 
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Appendix B. The disc load case 



(Jzzir) 



As an exercise and for verification purposes we describe a solution to the viscous half-space 
flat earth equation . The inverse Hankel transform of this solution is a solution to equation 

m- 

Consider a disc load centered at the map-plane origin, of radius Rq and thickness Hq. In 
this case 

-pigHo, <r < Rq, 
0, Ro <r 

where pi is the density of the load; we might as well suppose it is ice. Actually, let's suppose 
this load is applied at time zero and is held in place: azz{f-,t) = azz{r)H{t). Because this 
load is radial, the Hankel-transformed equation ((2)) is useful. We will assume an undeformed 
state vX = Q at time t = 0. 

We need the Hankel transform azz'- 

^zz{k) = - PigHo Jo{Kr)rdr = - ' ^ / Jo{s)s ds = -pigHoRon^^ Ji{kRo), 
Jo 1^ Jo 

using the change-of- variable s = nr and the identity ^(■sJi(s)) = sJo(s) (formula (8) in 
Appendix A of (Sneddon 1951)). 

Now, we can solve © because it is simply a collection of decoupled first-order linear ODE 
problems in time: 

{exp{-(3{n)t/i2r]K)) - 1} M^^Ro) 



u {K,t) = pigHoRo- 
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for t > and {K,t) = for t < 0, where P{k) = prg + Dk^. The displacement can be 
found by an (inverse) Hankel transform 

/•oo 

u^(r, t) = pigHoRo / Pi^)'^ {exp(-/?(K)t/(2r?K)) - 1} Ji{kRo) Mt^r) dn. (43) 
Jo 

This integral can be computed numerically, but care must be taken because the integrand 
is quite oscillatory. We break up the integrand into many (> 100) subintervals and call Mat- 
LAB's quadl on each subinterval. The code viscdisc in Appendix IdI implements equation 
(|43|1 . In particular, we believe that computing (|43|) for geophysically reasonable parameter 
values is quite accurate, and that the result can be used to verify the result from the Fourier 
spectral collocation method described in section |31 See section [5] for results. The displace- 
ment for a particular disc load is graphed in figure ITHl This solution looks rather like the 
Green's function plotted in figured as expected, but with a wider depressed area. 

The equilibrium limit of this disc load solution is of interest: 

/•oo 

u°°{r) = lim u^{r,t) = -pigHoRo / P{k)-^ Ji{hRo) M^r) dn. (44) 

This function satisfies the PDE Prgu°° + DA^u°° = a^z- It appears in figure El and we see 
the peripheral "dip" and "bulge" which occur at the edge of the disc. Note that the central 
part of disc load of thickness Hq and large radius descends to the "compensation depth" 
-{pi/pr)Ho. 
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Figure 16. Vertical displacement at 2000 years, and the equilibrium position 
at time oo, for a disc load of ice with thickness 1000 (m) and radius 1000 (km). 
"Compensation depth" corresponding to pi = 910, pr = 3300 (kgm^^) also 
shown. Note log scale on horizontal axis. 

Appendix C. On "Green's function thinking" for PDE (|T7jl 

Here we make some comments on formula © in section ^ which relate it to results in 
section IHl The time integral in © starts from — oo to avoid requiring precise knowledge 
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of the displacement {x,y,to) at any particular past time to- That is, © computes the 
displacement at time t caused by a load history known so far into the past that prior displace- 
ment states are irrelevant. This is reasonable because the underlying PDE, equation (|17|) in 
section 13 is diffusive and thus the influence of any displacement state decays exponentially 
in time. 

If, on the other hand, an initial displacement u{x, y, to) is known at a relatively recent past 
time then we must adjust As a start, integrating © by-parts in the time variable and 
assuming ^(x,y, — oo) is zero gives the formula 

^^(x,y,t)= /* ff^{\r-v%t-t')^{x',y\t')dx'dy'dt'. (45) 



— oo 



R 



Our assumption that there is no load at t = — oo is equivalent to assuming ^'(x,y, t) = 
A(x, y, x) ds. The significance of is that the load function ^ reappears and that we 
are motivated to examine dG^ /dt. 
From ((21) note that 

(r^t) = -^ {2r]Ky^exp[-(3{K)t/{2r]K)] Jo{Kr)KdK. (46) 







9t ' ' ' 27r 

In fact, formula (|46)) can also be extracted from equation (|2()j) in sectional 



dG^. , g ^_J exp[-/3(C,C)V(2r?(e^+0^/^)] l 

^(x,,,t) = --^, 2rjie + C^)V2 (4^) 



= / {2r,K)-^ exp [-(3{K)t/{2r]K)] J^i^x^ + y^ K)KdK. 

Here J-2 stands for the two-variable Fourier transform and I3{k) = prQ + Dk^ as in section 
n The second equality in (|47j) is explained by noting (Appendix E)) that becomes the 
Hankel transform on radial functions. Comparing to section ^ azz{x,y,t) = —g^{x,y,t) 
(N m-2) if $ (kg m ^) is the load function. 

Returning now to equation (|45|) . we can see from (|20|) how to modify (|45|) to include 
knowledge of displacement at a finite time to- In fact, the convolution theorem 

^2' {f9} = ^(f*9) = ^JJ fix -x',y- y')g{x', y') dx' dy' 

now allows us to write the inverse Fourier transform of (|2()j) as 

u^(x,y,t)= [' ll^{\r-r'\,t-s)^ix',y',s)dx'dy'ds (48) 



to 



+ // j{\r-r'\,t-to)uix',y',to)dx'dy' 



jir,t) = — / exp [— /3(K)t/(2?7K)] jQirK)KdK. 







where 

^(r. t) = 

2tx 

In our view equation (|48l) proves equation We have justified the heuristic Green's 
function thinking in section ^ by a standard linear analysis of PDE (|17() . Note that as 
t — > 00, 7(r, t) 0. Thus as to — > —00 the last term in (|l8|) disappears. The decaying 
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kernel 7 describes the rate at which information held in the previous displacement u{x, y, to) 
is eliminated. 

Appendix D. Matlab codes 

This section gives Matlab implementations of the numerical strategies outlined in sections 
121121 and AppendixEl Here is a synopsis of the three codes: 

• gef orconv precomputes the spherical elastic load response matrix given by equa- 
tion Q, which is the computational form of the Green's function . It numeri- 
cally computes integrals @ by using Matlab 's dblquad with the default settings, 
where the integrand comes from linearly interpolating the tabular data given by 
(Farrell 1972). If geforconv has been run on the given grid, fastearth (below) 
then convolves with the load, as noted, using Matlab's conv2 for now because it 
seems sufficiently fast; an FFT-based alternative is reasonable. The precomputation 
is much more expensive than running the viscous deformation model or convolving 
with the load. For instance, in the Matlab run 

tic, I = gef orconv (256, 256, 2000, 2000) ; toe 
save 1256 I 

tic, fastearth(256, 100000,500,4, 'teste) ; toe 



the call to geforconv requires about 4 hours while only 10 minutes is required to 
run fastearth; the machine is a 3 GHz Pentium IV. 

• fastearth computes the viscous model using f ft2 to compute equation H26|). It also 
convolves with the load for the spherical elastic model as noted. In this example 
implementation the load is either a standard disc load or a modification to Test C 
from (Bueler et al. 2005). 

• viscdisc implements formula (|43|) in Appendix IbI 

gef orconv .m: 



function II=gef orconv (Nx,Ny,Lx,Ly ) ; 

7, GEFORCONV Computes matrix I=I(p,q) , a form of the elastic "load response 
7, matrix" described in Lingle & Clark (1985) and in the technical report 
'/, Bueler (2005). Matrix l(p,q) is used by FASTEARTH. 
7. 

7, 1=GEF0RCQNV(NX,NY,LX,LY) computes (2 NX + 1) by (2 NY + 1) 

7o matrix I with entries 
7 /dy/2 /dx/2 

7 l(p,q)= I I G"E(sqrt((pdx-xi)"2+(qdy+eta"2))) dxi deta 

7 /-dy/2 /-dx/2 

7 where dx=2 LX/NX, dy=2 LY/NY. Input LX,LY in km. Assumes region 

7 extends from x=-LX to x=LX and y=-LY to y-LY. Uses Matlab's 
7 dblquad to do integral. 

7 

7 Example : 

7 l=geforconv(16, 16,2000,2000) ; 

7 

7 See also FASTEARTH. 
7 ELB 12/11/05. 

rm=[ 0.011 0.111 1.112 2.224 3.336 4.448 6.672 8.896 11 . 12 17 . 79 . . . 
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22.24 27.80 33.36 44.48 55.60 66.72 88.96 111.2 133.4 
222.4 278.0 333.6 444.8 556.0 667.2 778.4 889.6 1001.0 
1334.0 1779.0 2224.0 2780.0 3336.0 4448.0 5560.0 6672.0 7784.0 
10008.0] * le3; 7. converted to meters 
% GE /(10~12 rm) is vertical displacement in meters 
GE=[-33.64 -33.56 -32.75 -31.86 -30.98 -30.12 -28.44 -26.87 -25.41 
-21.80 -20.02 -18.36 -17.18 -15.71 -14.91 -14.41 -13.69 -13.01 
-12.31 -10.95 -9.757 -8.519 -7.533 -6.131 -5.237 -4.660 -4.272 
-3.999 -3.798 -3.640 -3.392 -2.999 -2.619 -2.103 -1.530 -0.292 
0.848 1.676 2.083 2.057 1.643]; 
% linearly extrapolate GE to r=0; GE(0) := -33.6488 
GE=[interpl(rm(l:2) ,GE(1:2) ,0.0, 'linear' , 'extrap') GE] ; 
rm=[0.0 rm] ; '/. length (rm)=length(GE) =42 

dx=2*Lx*1000/Nx; dy=2*Ly*1000/Ny ; '/. dimensions of load element 
II=zeros(2*Nx-l,2*Ny-l) ; 
% compute entries of II by dblquad, using quad method and default T0L=le-6 
for p=-Nx+l:Nx-l 

for q=-Ny+l:Ny-l 

II (p+Nx , q+Ny) =dblquad(aintegrand , -dx/2 , dx/2 , -dy/2 , dy/2) ; 

end 

end 

'/, nested function which is integrand in I(p,q) 
function z=integrand(xi,eta) 

r=sqrt((p*dx-xi) . "2+(q*dy-eta) ."2) ; 
z=zeros (length(r) ) ; 
for j j=l : length(r) 

if r(jj)==0.0 '/, treat normalization r as 11 m 

z(jj)=GE(l)/(rm(2)*lel2) ; 
elseif r(jj)>=rm(end) , z(jj)=0; 

else y, linearly interpolate to get GE value; then normalize 
ii=find(rm>r(jj) ,1) ; 2 <= ii <= 42 

z(jj)=interpl(rm(ii-l:ii) ,GE(ii-l:ii) ,r(j j))/(r(jj)*lel2) ; 

end 

end 

end 

end 



177.9 ... 
1112.0 ... 
8896.0 ... 



f as t earth. m: 



function [uv , ue , H , xs , ys] =f astearth (Nin , tf , dtyear , Z , type , graphf lag) 

■/, FASTEARTH Use the elastic-viscous earth model from 

7, Bueler (2006) "Computation of a viscoelastic deformable earth 

7. model for ice sheet simulation." 

7i to simulate the deformation of the earth under a disc load or under an 

y, artificial ice sheet load history. Uses a modified "load response 

7, matrix" approach for the elastic component. Uses a Fourier spectral 

y, method (Trefethen 2000) for the flat-earth, elastic plate over viscous 

y, half -space model. The PDE for this latter model is 

y, d/dt(2 eta Lap*{l/2} u) + rho g u + D Lap"2 u = sigma.zz. 

y. Time discretization of this PDE is Crank-Nicolson. The source of the 

% model is 

% Lingle fe Clark (1985) "A numerical model of interactions between a 
'/, marine ice sheet and the solid earth: Application to a West Antarctic 
y, ice stream," J. Geophysical Research 90 (CI) 1100—1114. 



COMPUTATION OF A VISCOELASTIC DEFORMABLE EARTH MODEL 



•/. [UV,UE,HH,XX,YY]=FASTEARTH(N,TF,DTYEAR,Z) computes and plots 

■/, deformation of the earth's surface under an ice disc load of radius 
■/, lOOOkm and thickness 1000m. Uses N grid points in each direction on 

Z a square 4000 km by 4000 km region R of physical interest. N must be 
% even; best choices are powers of 2. Z is a small whole number 

(Z=2 recommended) which determines the size of the computational 
'/, region Omega for the flat viscous model. Namely, Omega is a square 
■/, with side 4000*Z km, discretized by N*Z points in each direction. A 
'/, periodic boundary condition is applied at the edge of Omega. TF is 

# of years of the run with timestep DTYEAR. Elastic deformation is 
■/o performed on R and is only done if a file "I{N}.mat" can be found; 
■/. this can be generated using GEFORCONV; see examples. FASTEARTH 
■/, returns the gridded flat viscous deformation UV, the gridded elastic 

deformation UE, the thickness HH — all at the final time and the 

■/. grid itself in XX, YY (as from MESHGRID) ; all outputs in meters. 
% 

t [UV, UE.HH, XX, YY]=FASTEARTH(N,TF, DTYEAR, Z, 'disc') Same as above. 

•/. 

■/. [UV, UE,HH, XX, YY]=FASTEARTH(N,TF, DTYEAR, Z, 'teste) is the same except 

7, that the ice load is essentially test C, but with simple isostacy 
7, added, from 

7, Bueler, et al. (2005) "Exact solutions and the verification of 

% numerical models for isothermal ice sheets," J. Glaciol. 51 (173) 

■/. 291--306. 

7. Thickness evolution is stopped at 40033 years and the last thickness 

7. held. 

7. 

7. [UV, UE.HH, XX, YY]=FASTEARTH(N, TF, DTYEAR,Z, '...' ,false) No graphing. 

7. 

7, Example I; generate I for 16 x 16 case, store and use: 

7. tic, 1 = geforconv(16, 16, 2000, 2000) ; toe 7. about a minute 

7. save '116' I 

7. fastearth(16,100000,500,2) ; 7. 2 sees 

7i Example II; to see ice cap on bed; needs "1128. mat": 

7, [uv, ue,H,x,y]=fastearth(128, 40000, 500, 2, 'teste) ; 7. 15 sees 

7. figure, mesh(x/1000,y/1000,H+uv+ue) 

7. hold on, mesh(x/1000,y/1000,uv+ue) , hold off 

7. See also GEFORCONV, VISCDISC, TESTFAST. 

7. ELB 1/12/06 



if nargin<5, Cflag=false; type='disc'; end 
if nargin<6, graphf lag=true; end 
if stremp(type, 'teste) , Cflag=true; 
else, Cflag=false; end 

if ( (floor (Z)~=eeil(Z)) I I (Z<1)), errorCZ must be positive integer'), end 



g=9.81; spera=3 . 1556926e7; 7. seconds per year 

rhoi=0 . 910e3 ; 7. density of ice; kg/m~3 

rhor=3 . 300e3 ; 7 density of mantle; kg/m~3 

D=5.0e24; 7. flexural rigidity of lithosphere as an elastic plate 

eta=1.0e21; 7. viscosity of mantle 



L0=2000e3; 7. half-length of actual domain in each direction; m 

7o compute flat-earth, viscous half -space model on extended grid: 

L=Z*2000e3; 7. half-length of computational domain 

N=Z*Nin; 7. computational N 
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h=2*L/N; x=-L+h:li:L; 

[xx,yy]=meshgrid(x,x) ; rr=sqrt(xx. "2+yy. "2) ; 
M=ceil(tf/dtyear) ; dt=(tf/M)*spera; 

'/, Fourier coefficients of powers of Laplacian 
cx= (pi/L) * [0 : N/2 N/2-1 : -1 : 1] ; 
[ccxx,ccyy]=meshgrid(cx,cx) ; 
cclap=ccxx. "2+ccyy . "2; 

cchalf =sqrt (cclap) ; '/, lap"{l/2]- coeffs 
ccbih=cclap. "2; % biharmonic coeffs 

■/, coefficients in transformed PDE 

partl=2*eta*cchalf ; 

part2=(dt/2)*( rlior*g*ones(size(ccbih)) + D*ccbih ); 

right = parti - part2; % Fourier-transformed operator on right of eqn (19) 
left = parti + part2; % ... on left of eqn 

dispC computing flat viscous deformation ...') 
uunl=zeros (size (xx) ) ; 
for n=0:M-l 
uun=uunl ; 

if Cflag, H=getTestC(xx,yy, (n+l/2)*dt) ; 
else, H=1000*(rr<le6) ; end 
y, apply Fourier spectral method: 
sszz=-rhoi*g*H ; 

frhs=right.*fft2(uun) + f f t2(dt*sszz) ; 
uunl=real(if ft2( f rhs . /lef t )); 

end 

'/, tweak: find average value along "distant" bdry of [-ZL , ZL] X [-ZL , ZL] , 
■/, remove it and add back result for equivalent disc load 
uunl=uuni-( sum(uunl (1 , : ) )+suiii(uunl ( : , 1) ) )/(2*N); 
if Cflag '/, assume equiv disc load has R0=1000km 

H0=h*h*sum(sum(H))/(pi*le6*2) ; trapezoid rule 
else, H0=1000; end 

uunl=uunl+viscdisc(H0,1000,tf ,L/1000) ; 

sh=(N/2)-(N/(2*Z))+l: (N/2) + (N/(2*Z)) ; '/. extract central/actual part 
xs=xx(sh,sh) ; ys=yy (sh, sh) ; uv=uunl(sh,sh) ; H=H(sh,sh); 

'/, get and use (w. convolution) elastic LRM if available; plot things 

f ilename= [ ' I ' num2str(Nin) '.mat']; 
doGE= (exist (filename) ==2) ; 

if doGE, dispC [elastic load response matrix FOUND]') 
else, dispC [elastic load response matrix NOT found]'), end 
if doGE 

disp([' computing spherical elastic deformation using I' ... 

num2str (Nin) ' . mat ...']) 
S=load (filename) ; 
ue=rhoi*conv2(H,S . I, ' same' ) ; 

if graphflag 

figure(l), clf, subplot (2, 1 ,2) , mesh(xs/1000,ys/1000,ue) ; 

xlabeK'x (km)'), ylabelCy (km)') 

zlabelCvertical displacement (m)') 

title (' spherical , self -gravitating elastic model') 

subplot (2, 1,1) , mesh(xs/1000,ys/1000,uv) ; 

end 

else 

dispC skipping spherical elastic deformation ...') 
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ue=zeros(size(uv)) ; '/, for return value only 

if graphfleig, figure(l), clf, mesh(xs/1000,ys/1000,uv) ; end 

end 

if graphflag 

xlabelCx (km)'), ylabelCy (km)') 
zlabeK 'vertical displacement (m)') 
title ('elastic plate over viscous half -space model') 

end 

function H=getTestC(x,y ,t) 

% compute exact H from test C (with simple isostacy added) in Bueler et al 
"Exact solutions ... isothermal ice sheets". Stops growth at time 

7, t0=40033 years. 

spera=3. 1556926e7; '/, seconds per year 
H0=3600; R0=750000; 

see "Exact solutions ..." eqn (9) : 
•/. t0=(2/9.0177e-13)*(7/((l-(910/3300))*4))-3*(R0-4/H0-7) is 40034 years 
tO=40034*spera; 

ts=t/tO; if ts>l, ts=1.0; end 
rscl= (ts" (-2) ) * (sqrt (x . -2+y . "2) /RO) ; 
H=HO*ts*max(0, 1-rscl . * (4/3) ).-(3/7); 



viscdisc .m: 



function u=viscdisc(HO,ROkm,tyrs ,rkm) 

■/, VISCDISC Solution of viscous half-space overlain by rigid lithosphere in 

case of disc load at origin. Computed by Hankel transform means for 
'/, comparison to Lingle ft Clark (1985) and to new numerical method in 
% technical report: Bueler (2005), "Computation of a viscoelastic 
% deformable earth model for ice sheet simulation." 
•/. 

■/. [R,U]=VISCDISC(HO,ROKM,TYRS,RKM) computes the displacement U at 
■/, time TYRS (years) at radii RKM (km) caused by a disc of radius 
% ROKM (km) and thickness HO (m) . RKM may be a vector; U will be a 
■/, vector of the same length. 
7. 

7. [R,U]=VISCDISC(HO,ROKM, 'INF' ,RKM) computes the equilibrium 
7i displacement . 

7. 

% Example I; for disc 1000 m thick with lOOOkm radius at t=10~2 , 10"4, inf 
7. years; plotted on range 0.01 km to 10000 km. Each evaluation of 
7i viscdisc requires more than a minute on my machine. 
7. 

7. r=10.-(-2: .01:4); 7. 601 r values 
7. ul=viscdisc(1000, 1000, 100, r) ; 
7. u2=viscdisc (1000, 1000, 2000, r) ; 
7. u3=viscdisc(1000, 1000, 'inf ' ,r) ; 

7. semilogx(r ,ul,r,u2,r ,u3) , xlabeK'r (km)'), ylabelCu (m)') 

7. u3(l)/1000, -910/3300 % compensation depth comparison 

7. 

y. Example II; continues the above. Plots the Green's function as in figure 
7. 5(a) in Lingle ft Clark (1985), but at fewer times. Takes a minute or so. 
7. 

7. tt=[20 200 500 2000 5000 20000 50000 100000]; 
7. R0=0.001; H0=l/(pi*910) ; 7. small disc with mass 1 kg 
7. for j=l:8, uu(: ,j)=viscdisc(HO, 0.001, tt(j) ,r) ; end 
7. figure, for j=l:8, semilogx(r ,uu( : , j ) ) , hold on, end 
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'/, hold off, xlabeK'r (km)'), ylabelC displacement (m)') 

■/, axis([r(l) r(end) -3.5e-15 0.5e-15]), grid on 

•/. 

•/. See also GHV, FASTEARTH, TESTFAST. 
•/. ELB (12/7/05) 

g=9.81; % m/s"2 
rhoi=910; '/. kg/m"3 
rhor=3300; '/. kg/m"3 
rg=rhor*g; '/. combined constant 

D=5.0e24; 7, N m; flexural rigidity of lithosphere 
eta=1.0e21; % Pa s; viscosity of mantle 
spera=3. 1556926e7; '/. seconds per year 

R0=R0km*1000; 
TOL=eps ; 

pts=[10.-(-3:-0.05:-10) l.Oe-14] ; 

if (isstr(tyrs)&festrcmp(tyrs, 'inf ')) 
for k=l : length (rkm) 
rk=rkm(k)*1000; 

result=quadl(Sequilgrand,pts(l) , 100.0*pts(l) ,TOL) ; •/.kap->infty tail 
for j=l: length (pts)-l 

result=result+quadl(@equilgrand,pts(j+l) ,pts(j) ,TOL) ; 

end 

u (k) =-rhoi*g*HO*RO*result ; 

end 

else 

t=tyrs*spera; 

for k=l : length (rkm) 

rk=rkm(k)*1000; 

result=quadl(Sintegrand,pts(l) ,iOO.O*pts(l) ,TOL) ; 7, kap->infty tail 
for j=l: length (pts)-l 

result=result+quadl(9integrand,pts(j+l) ,pts(j) ,TOL) ; 

end 

u(k)=rhoi*g*HO*RO*result ; 

end 

end 

function y=integrand(kap) 

7. integrand of inverse Hankel transform 

beta=rg + D*kap."4; 

expdif f =exp (-beta*t . / (2*eta*kap) ) -ones (size (kap) ) ; 
y=expdif f . *bessel j (1.0,kap*R0) .*besselj (0.0,kap*rk) ./beta; 
end 

function y=equilgrand(kap) 

7. integrand of inverse Hankel transform when t — >infty 
y=besselj (1.0,kap*R0) .*besselj (0.0,kap*rk) ./(rg + D*kap."4); 

end 



end 



